Land Cover 2010 to 2020

Land Cover 2010 to 2020#

From 2010 to 2020, the Denver MSA has maintained control of urban sprawl. While the population increased by 16%, from 2.54 million to 2.95 million (block-level population from the 2010 and 2020 decennial census datasets), the urban area has expanded by only 3.5%, from 1,200 to 1,240 square miles, according to data from MRLC.

We used land cover data from 2008 and 2019 as proxies of the land cover in 2010 and 2020, as these are the best available data from MRLC. Land cover data rasters with 4000*4000 sqft cell sizes. Our study identified 5 types of land cover types:

Urban.
Forests.
Farms.
Wetlands.
Water.

Show code cell content Hide code cell content

# Read land cover rasters
# The file paths are relative to the notebook source file's location in the repo's main branch
land_cover_2009 = rio.open("../data/land-cover/2009.tif")
land_cover_2019 = rio.open("../data/land-cover/2019.tif")

def raster_to_shape(raster, crs, to_point=False, value_column="value"):
    data_array = raster.read(1)
    transform = raster.transform
    shapes_data = []
    height, width = data_array.shape

    if to_point:
        for row in range(height):
            for col in range(width):
                value = data_array[row, col]
                if np.isnan(value):
                    continue
                x, y = transform * (col + 0.5, row + 0.5)
                # Create a polygon for the cell
                point = Point(x, y)
                # Add the polygon and its value to the list
                shapes_data.append({"geometry": point, value_column: value})
        return gpd.GeoDataFrame(shapes_data, crs=crs)

    for row in range(height):
        for col in range(width):
            value = data_array[row, col]
            if np.isnan(value):
                continue
            x_min, y_max = transform * (col, row)
            x_max, y_min = transform * (col + 1, row + 1)
            # Create a polygon for the cell
            polygon = Polygon(
                [(x_min, y_min), (x_min, y_max), (x_max, y_max), (x_max, y_min)]
            )
            # Add the polygon and its value to the list
            shapes_data.append({"geometry": polygon, value_column: value})
    return gpd.GeoDataFrame(shapes_data, crs=crs)

# Get a fishnet of info from the land cover rasters
fishnet = (
    raster_to_shape(
        land_cover_2009, crs, to_point=False, value_column="land_cover_2009"
    )
    .query("land_cover_2009 != 0")
    .copy()
)
points_2019 = (
    raster_to_shape(land_cover_2019, crs, to_point=True, value_column="land_cover_2019")
    .query("land_cover_2019 != 0")
    .copy()
)
# Join points to fishnet, no index right
fishnet = gpd.sjoin(fishnet, points_2019, how="left", predicate="contains").drop(
    columns=["index_right"]
)

# Classify land cover types
developed = [21, 22, 23, 24]
land_covers = {
    "developed": [21, 22, 23, 24],
    "forest": [41, 42, 43],
    "farm": [81, 82],
    "wetland": [90, 95],
    "water": [11],
}
for year in [2009, 2019]:
    # First update whether developed bool
    fishnet[f"developed_{year}"] = fishnet[f"land_cover_{year}"].isin(developed)

    # Then update land cover
    fishnet[f"land_cover_type_{year}"] = "other"
    for land_cover, codes in land_covers.items():
        fishnet.loc[
            fishnet[f"land_cover_{year}"].isin(codes), f"land_cover_type_{year}"
        ] = land_cover

    fishnet.drop(columns=[f"land_cover_{year}"], inplace=True)

The above plots show that the Denver MSA’s urban growth was modest during 2010-2020. Most new developments were sporadically located on the fringes of existing development. We can also see a patter that developments tend to be located near some key corridors to the southeast of the City of Denver. We will further explore these patterns in the next sections.